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Quantum Monte Carlo methods provide in principle an accurate treatment of the many-body 
problem of the ground and excited states of condensed systems. In practice, however, uncontrolled 
errors such as those arising from the fixed-node and pseudopotential approximations often limit 
the quality of results. We show that the accuracy of quantum Monte Carlo calculations is limited 
by using available pseudopotentials. In particular, it is necessary to include angular momentum 
channels in the pseudopotential for excited angular momentum states and to choose the local channel 
appropriately to obtain accurate results. Variational and diffusion Monte Carlo calculations for Zn, 
O, and Si atoms and ions demonstrate that these issues can affect total energies by up to several eV 
for common pseudopotentials. Adding higher-angular momentum channels into the pseudopotential 
description reduces such errors drastically without a significant increase in computational cost. 
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I. INTRODUCTION 

Computational electronic structure methods have been 
extremely useful in developing our understanding of the 
atomic and electronic structures of real materials. As 
methods have become more accurate and their implemen- 
tations increasingly efficient, simulation and calculation 
have taken some of the burden of finding and character- 
izing new materials off experimental work. 1 

Density functional theory (DFT), in particular, has 
been widely applied to many systems in recent decades. 
It is computationally efficient compared to other methods 
with similar accuracy, and robust, user- friendly software 
packages have made the method easy to apply. How- 
ever, its accuracy is still insufficient for some applica- 
tions, and the lack of a systematic way to improve its 
results or estimate its errors has hampered progress. In 
particular, electron correlation effects can be significant 
in many complex materials and are not captured accu- 
rately by many of the commonly-used density function- 
als. The development of functionals which accurately 
describe the electronic band gap, van der Waals interac- 
tions, and other electronic properties of materials is still 
an active area of researchPHSl 

These issues are overcome by methods which treat 
quantum many-body effects explicitly from the outset 
such as quantum Monte Carlo (QMC). QMC methods 
are among the most accurate many-body methods and 
can reliably and accurately predict ground-state expec- 
tation values for many systems and have often been used 
as a benchmark for DFT workP Among the quantum 
Monte Carlo methods, variational Monte Carlo (VMC), 
diffusion Monte Carlo (DMC), and auxiliary- field quan- 
tum Monte Carld-^ are the most mature in terms of ap- 
plicability to solid state systems. We treat only VMC 
and DMC in this work and refer to them collectively as 



QMC. 

As computers become faster and high-quality soft- 
ware packages for QMC such as CASINO® QMCPACKP, 
QWALKP, and CHAMP 11 mature, these calculations 
are becoming less challenging. It is therefore important 
to identify and propagate the best-practice procedures 
for performing these calculations as they become more 
routine. 

QMC and other correlated-electron methods, particu- 
larly for heavy elements, usually employ the pseudopo- 
tential approximation to reduce the computational cost. 
The common form is the non-local, norm-conserving 
pseudopotential 12 which includes different potentials for 
each angular- momentum component of the wavefunction. 

In this work, we determine the error in the energy due 
to an insufficient number of angular-momentum channels 
in the pseudopotential and discuss other sources of error 
in QMC calculations. We show that pseudopotentials 
which include channels to account for higher angular- 
momentum components of the wavefunction are neces- 
sary for performing accurate pseudopotential calculations 
in QMC. Such pseudopotentials are not the norm in the 
literature, and we suggest that this be corrected in order 
that QMC methods be suitable for routine applications 
to scientifically and technologically interesting systems. 



II. BACKGROUND 



The computational cost of all-electron QMC scales ap- 
proximately as Z 5 - 5 or Z 6 - 5 with respect to the atomic 
number!- 13 * 14 * This scaling makes the direct application of 
all-electron QMC to heavy atoms difficult. In practice, 
many properties of atoms are primarily due to the be- 
havior of and interactions between valence electrons, and 
so a pseudopotential approximation is commonly used to 
remove core electrons from the calculation and reduce 



the necessary computational effort. 

Modern pseudopotentials are non-local in the sense 
that they act differently on distinct angular-momentum 
components of the wavef unction. This is necessary to ac- 
curately capture the effects of the nucleus and core elec- 
trons on the valence electrons since the pseudopotential 
not only represents the effective electrostatic potential 
but also enforces orthogonality of the valence orbitals to 
the lower-energy states of the same angular momentum. 

Of course, there is no clear distinction between core 
and valence electrons in many-body methods as the elec- 
trons are indistinguishable particles. Thus, the applica- 
tion of pseudopotentials does neglect exchange and cor- 
relation interactions between the valence and removed 
core electrons as well as that between the core electrons 
themselves. These errors are not explicitly accounted for 
in the calculations. However, the core-core interactions 
largely cancel out when considering energy differences, 
and the core-valence interactions may be kept small by 
a choice of core size which leads to significant spatial 
separation between the core and valence electron densi- 
ties. These techniques can lead to results as accurate as 
all-electron calculations.^ Additionally, the use of core- 
polarization potentials can account for some core- valence 
correlations GStCD 

To introduce non-local pseudopotentials in QMC, the 
electron-ion potential of any given atom is divided into a 
local potential Vj oc which is applied to the whole wave- 
function and several corrections V n i which account for 
the difference between the local potential and those which 
should be seen by the individual angular-momentum 
components of a wavefunction: 
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The nonlocal operator V K . acts on a function /(i"j) by 

I ' J 47T 

l,m 

(2) 
where the angular integral in the operator projects the 
wavefunction onto spherical harmonics. Each angular- 
momentum component thus "feels" its own potential, 

T)S 

V K (r) , which is only a function of the electron- nuclear 

distance r and accounts for the difference between the 
desired Z-dependent potential and the local channel V, . 

The local potential, or local channel, V, (r), is by con- 
vention chosen to be the exact potential applied to one 
of the angular- momentum components, and so sum in 
Eq. 12] need not include this local component P 

The choice of local channel itself is arbitrary but is 
often chosen for convenience during the pseudopoten- 
tial design process. In particular, judicous choice of the 
local channel is often necessary to avoid the problem 
of ghost states which can arise due to the Kleinman- 
Bylander transformation!^ We will see that the same 



choice of local channel that is suitable for that transfor- 
mation may not be the best with regards to accuracy of 
QMC calculations.^ 

In an independent electron theory such as Hartree- 
Fock or DFT, atomic wavefunctions are composed of 
some number of the lowest-energy single-particle or- 
bitals. For example, in these frameworks, the elec- 
tronic configuration of an oxygen atom may be writ- 
ten as ls 2 2s 2 2p 4 . Notice that this wavefunction con- 
tains no angular-momentum components above 1 = 1. 
Thus, a non-local pseudopotential in the above form 
which acts on these single-atom wavefunctions need not 

include terms V , for I > 2 if it is to be used to calculate 
nu — 

ground-state atomic properties. 

The situation is not so simple in the case of solids 
and other extended systems where changes in the wave- 
functions due to bonding effectively introduce higher 
angular-momentum components. Indeed, in the case of 
systems such as Si and other second row elements, one 
may find wavefunctions with significant higher angular- 
momentum character. In this case, it may be necessary 
to use a pseudopotential with a <i-channel when study- 
ing these systems in DFT, especially in the high-pressure 
regime. However, these errors often cancel when consid- 
ering en ergy differences and are frequently neglected in 
practicePEl 

In QMC and other correlated-electron methods, excita- 
tions of the wavefunction into higher angular-momentum 
states arise immediately, i.e., even for atoms. In VMC, 
wavefunctions may be represented by the product of a 
Slater determinant of single particle orbitals and the so- 
called Jastrow factor. The Jastrow is a positive function 
of inter-particle distances, and its purpose is to directly 
account for many-body correlation effects. Naturally, the 
VMC wavefunction is then no longer entirely composed 
of the lowest-order spherical harmonics. The situation 
in DMC (as well as other correlated-electron methods) is 
analogous P 

Notice from Eq. (II]) that, in the absence of a pseu- 
dopotential channel to deal with the higher-angular mo- 
mentum components of the wavefunction, these compo- 
nents simply feci the local channel. This is incorrect and 
may be drastically so, especially in the case where the 
local channel was designed to enforce orthogonality to 
the lower-energy orbitals with a particular angular mo- 
mentum. This can lead to sizeable errors in total energy 
calculations. 

Now, this effect is not a particularly surprising one and 
certainly has been understood by some in the density 
functional theory community since the beginnings of the 
use of pseudopotentials in that field (see, for example, 
Ref. [35]). However, inclusion of so-called higher angular- 
momentum channels is not the normal practice in the 
development of potentials for use with QMC. 

There are a limited number of pseudopotentials 
available for use with QMC. In particular, the ap- 
plication of projector-augmented waves^l or ultra-soft 
pseudopotential^S techniques in QMC is currently not 



TABLE I. Choices of angular-momentum channels and local 
channels for the various pseudopotentials considered for oxy- 
gen, silicon and zinc. 



Standard 
Channels Local 



Augmented 
Channels Local 



O 

Si 
Zn 



s, p 

s, p 

s, p, d 



s or p 
s or p 
s or d 



s, p, d 

s, p, d 

S, P, d, f 



s or p 

s or p 

s, d, or / 



feasible since the DMC operators for the projectors and 
the augmentation charge are unknown. However a num- 
ber of semi-local pseudopotentials have been developed 
with QMC applications in mind. Greeff et al. devel- 
oped a carbon pscudopotential which included s- and p- 
channelsPS Ovcharenko et al. applied a similar method- 
ology to produce pseudopotentials for Be to Ne and Al 
to Ar with Imax = 1 • Burkatzki et al. present poten- 
tials for many of the main group elements^ and for the 
3d transition metals™ Their Si and Zn potentials have 
3 channels, and their O potential has 2. These authors 
all cite the rule of thumb that Zmax should be at least 
as high as the highest angular-momentum component in 
atomic core. Trail et al. developed a variety of pseudopo- 
tentials for all elements from H to Hg. These all have 
exactly 3 channels and are associated with the CASINO 
code which, until recently, only supported pseudopoten- 
tials with exactly 3 channelsPS 



III. METHODOLOGY 

We determine how the number of channels and the 
choice of local channel affects the energy for several atoms 
and ions. We compute the total energies and first and 
second ionization energies of the zinc, oxygen, and silicon 
atoms using several related pseudopotentials. These el- 
ements provide interesting test cases due to their varied 
electronic structures. Additionally, we are interested in 
the application of Q MC m ethods to bulk semiconductors 
such as Si and ZnOJMMJ 

The oxygen and silicon pseudopotentials are based on 
those by Driver et aZp21, and the zinc pseudopotential is 
based on that by Bennett and RappePS All three are 
generated using the Opium pseudopotential codepS Cut- 
off radii and basis functions for the construction of the 
pseudopotentials were chosen to minimize the difference 
in pseudopotential and all-electron valence energy levels 
calculated in DFT using the PBE exchange-correlation 
functional^ for several electronic configurations. 

Table [T] lists the angular-momentum channels and the 
choice of local channel for each of our pseudopotentials. 
For each element, we consider (i) pseudopotentials with 
the minimum number of channels (s and p for Si and O; 
s, p and d for Zn) and (ii) pseudopotentials that con- 
tain an additional angular-momentum channel (d for O 
and Si; / for Zn). We refer to the first set as standard 
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FIG. 1. (color online) Pseudopotentials for O, Si, and Zn. 



pseudopotentials and the second one as augmented pseu- 
dopotentials. For the local channel we consider the s 
and p channels for Si and O and the s, d or / channels 
for Zn. This results in a total of 13 pseudopotentials 
as listed in Table [I] Figure [T] shows the distance depen- 
dence of the angular momentum channels for the various 
pseudopotentials. For Si and Zn, we confirmed that the 
pseudopotentials accurately describe the lattice parame- 
ters of the ground state crystal structure and for O, we 
confirmed that the pseudopotential reproduces the dimer 
bond length at the DFT level. 

QMC calculations were performed using the CASINO 
codeP We implemented support for pseudopotentials 
with an arbitrary number of angular-momentum chan- 
nels in CASINO. Total energy calculations are per- 
formed on the 9 isolated ions with each of the appli- 
cable pseudopotentials. The VMC calculations used 
Slater-Jastrow variational wavefunctions with orbitals 
expressed in a blip basis. 36 The single-particle orbitals 
were generated using the PWSCF code^ and the PBE 
exchange-correlation functional!^ Plane-wave cutoffs of 
70 Ry for oxygen and silicon and 100 Ry for zinc 
were used to converge the total energies to 2 meV. 
The known magnetic states of the atoms and ions were 
used. The Jastrow factor is a non-negative function of 
inter-particle distances and includes two-body electron- 
electron and electron-nucleus and three-body electron- 
electron- nucleus terms as implemented in CASINO™ 
The backfiow transformationPS was not found to pro- 



vide any significant benefit in these cases. Parameters 
were added to the Jastrow factor of the trial wavefunc- 
tion gradually during its optimization. The Jastrow pa- 
rameters were optimized using variance minimizations^ 
followed by energy minimization in the final steppH Trial 
wavefunctions were evaluated by their mean energy plus 
three times the statistical error in the energy, following 
Ref.Sl 

Several additional details of our VMC calculations are 
noteworthy. First, the integral in Equation pi) is per- 
formed on a spherical grid in real space. This integra- 
tion mesh must be chosen to be sufficiently dense to ac- 
curately calculate the contributions to the energy from 
higher angular- momentum compontents of the wavefunc- 
tion and thus evaluating the effects which are the focus 
of this paper. Secondly, it is the default behavior of the 
CASINO code that the non-local contributions to the 
energy are assumed constant and are not re-calculated 
during a variance minimization step. In many systems, 
this improves the runtime of the algorithm significantly 
while still giving good results — in some cases it actually 
improves the performance of the variance minimization. 
However, as we will see, the non-local contributions are 
significant in many of our calculations. We found it nec- 
essary in many cases to re-calculate the non-local con- 
tributions to the energy at each step of the optimization 
to ensure the stability of the optimization process during 
Jastrow optimization. 

Our DMC calculations were performed using the pseu- 
dopotential locality approximation.^ For each system, 
we performed 1,000 equilibration and 3,000 accumula- 
tion steps on each of 256 processors with a timestep of 
0.01 Ha -1 and a target population of 2,000 walkers. 

Finally, atomic ionization energies are simply dif- 
ferences between the total energies of the appropriate 
species. 



IV. RESULTS AND DISCUSSION 
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FIG. 2. (color online) Total VMC energy in Hartree and 
statistical error in the energy of each species with respect to 
each Hamiltonian. Pseudopotentials are denoted according to 
the choice of local channel and as 'aug' if they are augmented 
with an additional channel or 'std' otherwise. 
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Figures [2] and [3] show the energies for each ion- 
pseudopotential combination for VMC and DMC, respec- 
tively. The error bars indicate only the statistical uncer- 
tainties in the energies associated with the QMC calcu- 
lations. First, it is important to notice the axes. The 
variation in the total energies differs between the species. 
For Si, it is on the order of milli-Hartrees, while for Zn, it 
is on the order of tenths of Hartrees. O falls somewhere 
in between. 

Consider now the VMC and DMC total energies. The 
two sets of data exhibit similar trends which is to be ex- 
pected since they arise from the differences in the pseu- 
dopotentials. The first thing to notice is that the calcula- 
tions using an augmented pseudopotential (i.e., the two 
or three left-most data points in each panel) are largely in 
agreement with each other, while this is not the case for 
the standard potentials. That is, the choice of local chan- 
nel has a large effect on the total energy when using the 



FIG. 3. (color online) Total DMC energy in Hartree and 
statistical error in the energy of each species with respect to 
each Hamiltonian. Pseudopotentials are labeled as in Fig. [2] 



standard pseudopotentials since the higher-/ components 
of the wavefunction see that local channel. When using 
the augmented potentials, more of the wavefunction sees 
its correct potential, and the choice of local channel has 
less effect on the result of the calculation. 

Indeed, if we take the calculations with augmented po- 
tentials to indicate the correct result, we can understand 
the errors in the other total energies in terms of which 
potentials are incorrectly applied to certain components 
of the wavefunction. As seen in Figure [Tl the s-channel is 
the most repulsive for each of the species, the p-channel 
is in the middle, and the d-channel is the most attractive. 
The /-channel is slightly above <i-channel in the case of 
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TABLE II. Lowest-energy excitations in eV to higher-Z states 
for each species from experiment ! 44 " 47 l 
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FIG. 4. Variance in local energy in atomic units and associ- 
ated statistical error of each VMC calculation. Hamiltonians 
are labeled as in Figure [2] 



1st Ionization Energies 2nd Ionization Energies 







FIG. 5. Comparison of the ionization energies in eV for oxy- 
gen, silicon and zinc in DFT, VMC and DMC for the differ- 
ent choices of pseudopotential with experiment. Energies for 
oxygen and silicon include corrections for the pseudopotential 
error at the DFT level as described in the text. 



zinc. 

Thus, we expect that calculations in which components 
of the wavefunction incorrectly see the s-channel to be 
too high in energy. Indeed these data points (which are 
the second right-most point in each frame of the total en- 
ergy plots) exhibit this trend. Similarly, the right-most 
data point in each frame corresponds to a calculation 
wherein any higher I component of the wavefunction sees 
the (i-channel, and these results are erroneously low in 
energy. Even the residual differences between the ener- 
gies calculated using the augmented potentials follow this 
trend. This is indicative of small amounts of yet-higher 



Species 



O 



Si 



Zn 



Neutral 12.08 5.86 8.53 

Singly-Ionized 28.7 9.84 14.54 

Doubly-Ionized 40.23 17.72 31.9 



I character in the wavef unctions. 

We can get some idea of the amount of virtual exci- 
tations that might be present in the many-body ground 
state of each of the species by considering the energies 
of these excitations. In this way, we can understand the 
expected magnitude of these effects. 

Table |TT] presents the lowest energy excitation to a 
higher angular-momentum state for each of the three 
atoms for the various charge states. The excitation en- 
ergies of these states increase with the level of ioniza- 
tion. Additionally, the d levels are relatively high in 
oxygen but low in silicon. The / levels in Zn trend in 
between. Thus, we expect that the effects in total energy 
described in this paper will be especially pronounced for 
the neutral species relative to the positive ones and for 
silicon relative to oxygen. Note that this effect due to the 
lower excitation energies here has implications not only 
for the atomic wavefunctions. Lower-energy states are 
more likely to participate in bonding in molecules and 
solids, and it is especially important to design pseudopo- 
tentials to account for that. 

Indeed, for the silicon species, the decreasing signifi- 
cance of the extra channel with increasing charge is clear. 
This effect is less readily apparent in oxygen and zinc 
data, and is likely obscured by another notable and cor- 
related effect. The variance in the local energy for each 
VMC calculation is shown in Figure [4] Eigenstates of 
the Hamiltonian have a variance of zero, and higher vari- 
ances indicate worse approximations of the ground-state 
wavefunction. That is, higher variances are correlated 
with higher total energies. From the figure, it is clear 
that our calculations using the standard pseudopotentials 
with the s-channel local resulted in poorer-quality wave- 
functions. These higher variances are not fundamental 
properties of the system but demonstrate that it is more 
difficult to optimize a Slater- Jastrow wavefunction with 
respect to these less realistic Hamiltonians. Additionally, 
the wavefunctions might be better described by a multi- 
determinant expansion, especially in the case of oxygen. 

Finally, the first and second ionization energies for each 
element are shown in Figure [5j The DFT ionization en- 
ergies were computed from the same calculations used 
to create the orbitals for the VMC trial wavefunction, 
and finite-size effects due to periodic boundary conditions 
were treated using the method of Makov and Payne™ 
Several sources of errors may explain the deviation of 
these results from experimental values. First, we calcu- 
lated spin-orbit corrections to the total energies at the 
DFT level and found that they largely cancel in the ion- 



ization energies, resulting in corrections too small to ac- 
count for the observed differences in the ionization ener- 
gies between QMC, DFT and experiment. 

Second, the pseudopotential approximation itself leads 
to several errors other than those focused on in this pa- 
per. By removing explicit treatment of core electrons 
from the calculation, we are neglecting correlations be- 
tween the core and valence electrons. This is minimized 
but not altogether eliminated by designing pseudopoten- 
tials so that the core and valence electrons are spatially 
separated. The core-valence correlation may be particu- 
larly important for the case of zinc where the 3d valence 
electron state have a sizeable spatial overlap with the 3p 
core electron states and may explain the large errors in 
the ionization energies. Thirdly, evaluation of the pseu- 
dopotentials in DMC is subject to the locality approxima- 
tion^ used in this work or the lattice-regularized method 
by CasulaP^ Pozzo and AUSP™ found that, in magnesium 
and magnesium hydride, the errors of the locality approx- 
imation and the lattice-regularized method are compara- 
bly small, but that the lattice method requires a much 
smaller DMC time step. 

Finally, our pseudopotentials themselves could likely 
be further optimized within the same framework. In par- 
ticular, the scattering properties of PBE pseudopoten- 
tials may be poor at certain energy scales, and HF poten- 
tials might perform better in conjunction with correlated- 
electron methods.^ To test the accuracy of the pseu- 
dopotential approximation for oxygen, we performed an 
all-electron, single-determinant DMC calculation of an 
isolated oxygen atom with a Slater- Jastrow trial wave- 
function and found an ionization energy of 13.611(7) eV, 
in close agreement with the experimental value. This 
strongly suggests that the errors in the DMC ionization 
energies in Fig. [5] are due to the pseudopotentials. 

With this in mind, we estimated corrections for the 
pseudopotential error at the DFT level for oxygen and 
silicon. All-electron DFT/PBE ionization energies were 
calculated using the Gaussian code^ converged with re- 
spect to the atomic basis. The difference between these 
ionization energies and those found with PBE using pseu- 
dopotentials should capture much of the error due to the 
pseudopotentials, and we have added these differences to 
the QMC results for oxygen and silicon shown in Figure 

El 

In the case of zinc, the error in the ionization energy in 
QMC stems from the poor description of the 3d-levels of 
the zinc atom in DFT. Semilocal functionals are known 
to place the 3d level of the Zn atom significantly too 
hig] j52 | 53 j This results in an incorrect description of the 
rf-channel of the pseudopotential and of the 3d-orbital in 
the trial wave function which is reflected in both the large 
energy variance and large deviation of the QMC ioniza- 



tion energy from experiment. Correcting for this error at 
the DFT level by adding the energy difference between 
an all-electron and pseudopotential DFT calculation to 
the QMC results is found to be insufficent. 
V. CONCLUSIONS 

We showed that pseudopotentials which include chan- 
nels to account for higher angular-momentum compo- 
nents of the wavefunction are necessary for performing 
accurate pseudopotential calculations in QMC. For O, 
Si and Zn, we determined how the number of angular- 
momentum channels and the choice of local channel in 
the pseudopotential affects the total energy and ioniza- 
tion energies of these atoms in QMC. We find a sizable 
error in the total energies for any choice of local channel 
when the pseudopotentials do not include at least one 
additional angular-momentum channel above the high- 
est angular-momentum component of the ground state 
wavefunction of the atom. This is because, contrary to 
single-electron mean- field methods such as DFT and HF, 
atomic ground state wavefunctions in correlated-electron 
methods include higher angular-momentum character. 
These components effectively see the wrong potentials 
when using standard pseudopotentials. This situation is 
expected to be even more pronounced in the case of solids 
and molecules. 

Our results suggest that the best practice is to in- 
clude at least one channel in the pseudopotential above 
the highest angular-momentum component of the ground 
state wavefunction in single-particle methods. Addition- 
ally, this highest channel should be used as the local 
channel as it will generally be most similar to missing, 
yet-higher angular-momentum channels. 
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